"""
模块：driver_risk_system
功能：商用车驾驶行为分析与UBI风险建模
核心业务：
  1. 基于时空数据识别高风险驾驶行为
  2. 构建动态风险评分模型
  3. 为保险公司提供风险分群依据
依赖库：pandas, numpy, scikit-learn, geopandas, matplotlib
"""
import pandas as pd
import numpy as np
from sklearn.cluster import DBSCAN  # 基于密度的聚类
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt
from datetime import datetime, timedelta
import geopandas as gpd  # 地理空间分析

# ==== 1. 模拟数据集（20辆车5天的行驶数据） ====
# 设计要点：模拟真实业务场景中的GPS轨迹、速度波动和驾驶事件
np.random.seed(42)
num_vehicles = 20
days = 5
records_per_day = 288  # 5分钟间隔 (24h * 60 / 5)
total_records = num_vehicles * days * records_per_day

data = {
    "vehicle_id": np.repeat([f"V{str(i).zfill(3)}" for i in range(1, num_vehicles + 1)], days * records_per_day),
    "timestamp": pd.date_range("2025-07-01 00:00:00", periods=total_records, freq="5min"),
    # 速度模拟：正常分布+随机波动（单位：km/h）
    "speed": np.clip(np.random.normal(65, 15, total_records) + np.random.uniform(-10, 10, total_records), 0, 120),
    # 加速度模拟：95%正常范围，5%急加速/急刹车（单位：m/s²）
    "acceleration": np.where(
        np.random.rand(total_records) > 0.05,
        np.random.normal(0, 0.8, total_records),
        np.random.choice([2.8, -3.2], total_records)
    ),
    # 地理位置模拟：生成连续轨迹
    "gps_lat": np.cumsum(np.random.normal(0, 0.001, total_records)),
    "gps_lon": np.cumsum(np.random.normal(0, 0.001, total_records)),
    "driver_id": np.random.choice([f"D{str(i).zfill(3)}" for i in range(1, 31)], total_records)
}
df = pd.DataFrame(data)

# ==== 2. 数据预处理 ====
# 目标：清洗异常数据，提取时空特征
# 规则参考：ISO 39001道路交通安全管理标准

# 2.1 异常值清洗
df = df[
    (df["acceleration"].abs() <= 3.5) &  # 加速度阈值（物理极限）
    (df["speed"] >= 10)  # 过滤停车数据
]

# 2.2 时间特征工程
df["hour"] = df["timestamp"].dt.hour
df["is_night"] = np.where((df["hour"] >= 22) | (df["hour"] <= 5), 1, 0)  # 夜间驾驶标记
df["is_rush_hour"] = np.where(df["hour"].isin([7, 8, 17, 18]), 1, 0)  # 早晚高峰标记

# 2.3 连续驾驶时长计算（疲劳驾驶依据）
df["drive_duration"] = df.groupby("driver_id")["timestamp"].diff().dt.total_seconds() / 3600  # 单位：小时

# ==== 3. 高风险行为识别 ====
# 业务规则来源：中国《道路运输车辆动态监督管理办法[](@replace=10001)》

# 3.1 急加速：加速度 > 2.5 m/s² 持续5秒以上
df["hard_accel_flag"] = np.where(df["acceleration"] > 2.5, 1, 0)
# 标记连续事件（窗口滚动检测）
df["hard_accel_5s"] = df.groupby("vehicle_id")["hard_accel_flag"].rolling(6, min_periods=1).sum().reset_index(level=0, drop=True)
df["hard_accel"] = np.where(df["hard_accel_5s"] >= 5, 1, 0)

# 3.2 急刹车：加速度 < -2.5 m/s² 或速度突降
df["hard_brake"] = np.where(
    (df["acceleration"] < -2.5) | 
    (df["speed"].diff().abs() > 20),  # 5分钟内速度变化>20km/h
    1, 0
)

# 3.3 疲劳驾驶：夜间连续驾驶 > 4小时
df["fatigue_drive"] = np.where(
    (df["drive_duration"].fillna(0) > 4) & (df["is_night"] == 1),
    1, 0
)

# ==== 4. 动态风险评分模型 ====
# 模型原理：加权事件频率 + 空间风险密度

# 4.1 按驾驶员聚合风险事件
driver_risk = df.groupby("driver_id").agg(
    total_hard_accel=("hard_accel", "sum"),
    total_hard_brake=("hard_brake", "sum"),
    total_fatigue=("fatigue_drive", "sum"),
    night_drive_ratio=("is_night", "mean")  # 夜间驾驶占比
).reset_index()

# 4.2 空间风险分析（高危区域识别）
# 使用GeoPandas计算驾驶点与高危区域（如学校、施工区）的距离
gdf = gpd.GeoDataFrame(
    df, 
    geometry=gpd.points_from_xy(df["gps_lon"], df["gps_lat"])
)
# 模拟高危区域（真实场景需GIS数据导入）
danger_zones = gpd.GeoSeries([
    Point(116.40, 39.90),  # 模拟北京某学校
    Point(116.41, 39.91)   # 模拟施工区域
])
# 计算最近高危区域距离
gdf["min_danger_dist"] = gdf.geometry.apply(
    lambda p: danger_zones.distance(p).min()
)
# 合并到驾驶员数据
driver_risk = driver_risk.merge(
    gdf.groupby("driver_id")["min_danger_dist"].mean().reset_index(),
    on="driver_id"
)

# 4.3 风险评分计算（归一化0-100分）
weights = np.array([0.3, 0.4, 0.2, 0.1])  # 权重：急刹车>急加速>疲劳>高危区域
risk_components = driver_risk[[
    "total_hard_accel", "total_hard_brake", 
    "total_fatigue", "min_danger_dist"
]]
# 标准化处理（消除量纲影响）
scaler = StandardScaler()
scaled_components = scaler.fit_transform(risk_components)
# 加权计算原始风险分
driver_risk["raw_risk_score"] = np.dot(scaled_components, weights)
# 归一化到0-100分
driver_risk["risk_score"] = (driver_risk["raw_risk_score"] - driver_risk["raw_risk_score"].min()) / \
                            (driver_risk["raw_risk_score"].max() - driver_risk["raw_risk_score"].min()) * 100

# ==== 5. 风险分群（DBSCAN聚类） ====
# 选择DBSCAN原因：自动识别异常点（极端高风险驾驶员）
cluster_data = driver_risk[["risk_score", "night_drive_ratio", "min_danger_dist"]]
scaler = StandardScaler()
cluster_scaled = scaler.fit_transform(cluster_data)

# 聚类建模（参数说明：eps=邻域半径，min_samples=最小簇样本数）
dbscan = DBSCAN(eps=0.8, min_samples=3)
driver_risk["risk_cluster"] = dbscan.fit_predict(cluster_scaled)
# 映射业务标签
cluster_map = {
    -1: "极端风险",
    0: "低风险",
    1: "中风险",
    2: "高风险"
}
driver_risk["risk_level"] = driver_risk["risk_cluster"].map(cluster_map)

# ==== 6. 可视化与业务输出 ====
plt.figure(figsize=(15, 6))

# 6.1 风险等级分布（子图1）
plt.subplot(131)
risk_counts = driver_risk["risk_level"].value_counts()
colors = ["green", "blue", "orange", "red"]
plt.bar(risk_counts.index, risk_counts.values, color=colors)
plt.title("驾驶员风险等级分布", fontsize=12)
plt.xlabel("风险等级")
plt.ylabel("驾驶员数量")

# 6.2 风险评分箱线图（子图2）
plt.subplot(132)
driver_risk.boxplot(column="risk_score", by="risk_level", grid=False, vert=False)
plt.title("风险等级评分分布", fontsize=12)
plt.xlabel("风险评分")
plt.ylabel("")
plt.suptitle("")  # 清除自动标题

# 6.3 地理风险热力图（子图3）
plt.subplot(133)
high_risk_points = gdf[gdf["driver_id"].isin(
    driver_risk[driver_risk["risk_level"] == "极端风险"]["driver_id"]
)]
plt.scatter(
    gdf["gps_lon"], gdf["gps_lat"],
    c="lightgray", alpha=0.1, label="普通轨迹"
)
plt.scatter(
    high_risk_points["gps_lon"], high_risk_points["gps_lat"],
    c="red", alpha=0.5, label="高风险点"
)
plt.scatter(
    danger_zones.x, danger_zones.y,
    marker="X", s=100, c="black", label="高危区域"
)
plt.title("高风险驾驶空间分布", fontsize=12)
plt.legend()

plt.tight_layout()
plt.savefig("driver_risk_analysis.png")

# ==== 7. 业务策略生成 ====
# 7.1 实时告警规则（推送至车联网平台）
alert_rules = {
    "hard_accel": "急加速事件连续触发5次 → 推送‘平稳起步’提示",
    "hard_brake": "急刹车事件+速度突降 → 推送‘预判路况’提示",
    "fatigue_drive": "夜间连续驾驶>4小时 → 强制休息通知",
    "danger_zone": "进入高危区域500米范围 → 激活限速模式"
}

# 7.2 UBI保险定价策略（输出至保险公司系统）
ubi_strategy = driver_risk.groupby("risk_level").agg(
    avg_risk_score=("risk_score", "mean"),
    driver_count=("driver_id", "count")
).reset_index()
ubi_strategy["premium_adjust"] = ubi_strategy["avg_risk_score"].apply(
    lambda x: 0.85 if x < 30 else 1.0 if x < 70 else 1.3
)